#include <math.h>
#include "eno.h"
#include "fluids.h"

#include "eno.xtn"
#include "fluids.xtn"

extern void invert_matrix();

/* compute the stream function given the vorticity.
   (solve poisson problem)
   output in pcg_result.
*/
void init_stream_function()
{
int ix,iy;
int iterations;
double dz,dr,tolerance,coeff1;

  for_all_points
   if (point1_ok)
    {
     a00[ix][iy]=a01[ix][iy]=a02[ix][iy]=0.0;
     a10[ix][iy]=a11[ix][iy]=a12[ix][iy]=0.0;
     a20[ix][iy]=a21[ix][iy]=a22[ix][iy]=0.0;

     coeff1=-1/(2.0*sqr(h));
     a00[ix][iy]=coeff1;
     a20[ix][iy]=coeff1;
     a02[ix][iy]=coeff1;
     a22[ix][iy]=coeff1;

     a11[ix][iy]=-4.0*coeff1;

#define force_u(incx,incy) (gru(ix+incx,iy+incy))
#define force_v(incx,incy) (grv(ix+incx,iy+incy))


     dz=(force_u(0,1)+force_u(1,1)-force_u(0,0)-force_u(1,0))/(2.0*h);
     dr=(force_v(1,0)+force_v(1,1)-force_v(0,0)-force_v(0,1))/(2.0*h);

     bij[ix][iy]=(dr-dz);
    
     pcg_result[ix][iy]=s[ix][iy];  /* initial guess */ 
    }    /* point1_ok - looping all points on the STAGGERED GRID */

  tolerance=1e-7;

  printf("p-stream\n");
  top_stream_value_global=top_stream_value;
  bottom_stream_value_global=bottom_stream_value;
  invert_matrix(tolerance,&iterations,pcg_result);
  top_stream_value_global=bottom_stream_value_global=0.0;
  printf("it %4d",iterations);
}   /* init_stream_function */
